///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file restricted_pairing_function.h
///
/// \brief   header file for the restricted pairing function class
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef PCPS_RESTRICTED_PAIRING_FUNCTION_HEADER
#define PCPS_RESTRICTED_PAIRING_FUNCTION_HEADER

#include<src/pcps.h>
#include<src/input.h>
#include<src/coeff_obj.h>
#include<src/mantissa_exponent.h>
#include<src/wfn_data.h>
#include<src/command.h>

namespace pcps {

  template <class S> class Hamiltonian;

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   Class to hold a description of the restricted pairing function for a specific
  ///          configuration of the lattice.
  ///          
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template<class S> class RPF_Descriptor {

    private:

      /// \brief  the number of alpha sites
      int _nas;

      /// \brief  the number of alpha particles
      int _nap;

      /// \brief  the newly occupied site upon updating to the next configuration in the heirarchy
      int _new_occ_site;

      /// \brief  the newly empty site upon updating to the next configuration in the heirarchy
      int _old_occ_site;

      /// \brief  the list of occupied alpha sites ordered to match the rows of the occupied pairing matrix
      std::vector<int> _a_sites;

      /// \brief  the list of occupied beta sites ordered to match the columns of the occupied pairing matrix
      std::vector<int> _b_sites;

      /// \brief  the row numbers of the occupied pairing matrix, ordered by their corresponding sites' indices
      std::vector<int> _a_ordered_perm;

      /// \brief  the column numbers of the occupied pairing matrix, ordered by their corresponding sites' indices
      std::vector<int> _b_ordered_perm;

      /// \brief  for each row of the occupied pairing matrix, the corresponding position in _a_ordered_perm
      std::vector<int> _a_perm_pos;

      /// \brief  for each column of the occupied pairing matrix, the corresponding position in _b_ordered_perm
      std::vector<int> _b_perm_pos;

    public:

      int new_occ_site() const { return _new_occ_site; }
      int old_occ_site() const { return _old_occ_site; }
      int & new_occ_site() { return _new_occ_site; }
      int & old_occ_site() { return _old_occ_site; }

      RPF_Descriptor(const pcps::Input & userinp);

      void initialize_permutation(const std::vector<int> & sites, std::vector<int> & ordered, std::vector<int> & perm_pos);
      void prepare_occ_info(const pcps::Configuration<S> & config);
      void update(int new_ind, int old_ind);
      void build_occupied_matrix(const std::vector<S> & pm, std::vector<S> & opm) const;
      int get_row_col_index(int site) const;
      void get_row_col_diff(const std::vector<S> & pm, S * const vec,
                            int new_ind = -1, int old_ind = -1) const;
      double sign_change() const;
      double reordering_sign() const;
      void printout(const std::string & name) const;

  };

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   Class to represent a restricted pairing function.  Here restricted means that each
  ///          pair of spins must have an alpha and a beta, not that the overall wavefunction is
  ///          necessarily a singlet.  If the pairing matrix is symmetric, however, then the
  ///          the wavefunction is indeed a good singlet.
  ///
  ///          The pairing matrix rows correspond the alpha particles, while the columns correspond
  ///          to beta particles.
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template<class S> class RestPairingFunc : public pcps::CoeffObj<S> {

    private:

      /// \brief  the number of alpha sites
      int _nas;

      /// \brief  the number of alpha particles
      int _nap;

      /// \brief  the number of unoccupied alpha sites
      int _nau;

      /// \brief  the site to remove a particle from to get to the best resting configuration to pair with the
      ///         proposed new active configuration
      int _best_prop_p;

      /// \brief  the site to add a particle to to get to the best resting configuration to pair with the
      ///         proposed new active configuration
      int _best_prop_q;

      /// \brief  whether we are sampling w.r.t. to the wavefunction or the local rdm
      bool _wfn_sampling;

      /// \brief  whether to weight the local rdm elements by the [pq,pq] integrals
      bool _pqpq_guiding;

      /// \brief  how much to damp the diagonal contributions to the rdm importance sampling function
      double _diag_rdm_damping;

      /// \brief  ratio between the determinant value at the best resting configuration to pair with the
      ///         proposed new active configuration and the current resting configuration
      S _best_prop_ratio;

      /// \brief  the value of the rdm guiding function at the active configuration
      S _acti_ratio_sum;

      /// \brief  the value of the rdm guiding function at the most recently proposed configuration
      S _prop_ratio_sum;

      /// \brief  the value of the pairing matrix determinant at the most recently initialized configuration
      pcps::MantExp<S> _init_value;

      /// \brief  the value of the pairing matrix determinant at the resting configuration
      pcps::MantExp<S> _rest_value;

      /// \brief  the alpha first excitation ratios at the active configuration
      std::vector<S> _acti_a_ratios;

      /// \brief  the beta first excitation ratios at the active configuration
      std::vector<S> _acti_b_ratios;

      /// \brief  the alpha first excitation ratios at the most recently proposed configuration
      std::vector<S> _prop_a_ratios;

      /// \brief  the beta first excitation ratios at the most recently proposed configuration
      std::vector<S> _prop_b_ratios;

      /// \brief  the inverse of the occupied pairing matrix at the resting configuration
      std::vector<S> _Zmat;

      /// \brief  the transpose of the inverse of the occupied pairing matrix at the resting configuration
      std::vector<S> _ZTmat;

      /// \brief  the U matrix used in the matrix determinant lemma and woodbury matrix identity
      std::vector<S> _Umat;

      /// \brief  the V matrix used in the matrix determinant lemma and woodbury matrix identity
      std::vector<S> _Vmat;

      /// \brief  the VZ matrix used in the matrix determinant lemma and woodbury matrix identity
      std::vector<S> _VZmat;

      /// \brief  the UZT matrix used in the matrix determinant lemma and woodbury matrix identity
      std::vector<S> _UZTmat;

      /// \brief  the Q matrix used in the matrix determinant lemma and woodbury matrix identity
      std::vector<S> _Qmat;

      /// \brief  the identity matrix used in the matrix determinant lemma and woodbury matrix identity
      std::vector<S> _Imat;

      /// \brief  floating point workspace
      std::vector<S> _work;

      /// \brief  [pq,pq] two electron integrals
      std::vector<S> _pqpq_tei;

      /// \brief  the Umask array used in the matrix determinant lemma and woodbury matrix identity
      std::vector<int> _Umask;

      /// \brief  the Vmask array used in the matrix determinant lemma and woodbury matrix identity
      std::vector<int> _Vmask;

      /// \brief  the Qmask array used in the matrix determinant lemma and woodbury matrix identity
      std::vector<int> _Qmask;

      /// \brief  integer workspace
      std::vector<int> _iwork;

      /// \brief  description of the pairing function at the resting configuration
      pcps::RPF_Descriptor<S> _rest_desc;

      /// \brief  description of the pairing function at the active configuration
      pcps::RPF_Descriptor<S> _acti_desc;

      /// \brief  description of the pairing function at the proposed configuration
      pcps::RPF_Descriptor<S> _prop_desc;

      /// \brief  map used for looking up correct values for n2 in sto3g for testing
      std::map<std::string, std::pair<double, double> > _cmap;

      //void compute_noise_products(const pcps::Configuration<S> & config);
      //void build_little_mat(const pcps::ConfigDiff<S> & diff);
      //pcps::MantExp<S> get_single_excite_value(int new_ind, int old_ind, const std::vector<S> & m) const;
      //void print_occupancy_info();
      //void swap_row_col(const int b, const int i, const int j);
      //void test_noise_removal(const pcps::Configuration<S> & config);
      //void get_single_excitation_a_ratios(const pcps::Configuration<S> & config, const int b, std::vector<S> & a_ratios);

      void get_update_vec_V_U(const int vp, const int b, const pcps::RPF_Descriptor<S> & d, const int rc_ind);
      void get_update_vec_UZT_VZ(const int vp, const int b, const int rc_ind);
      void get_update_vecs(const int vp, const pcps::RPF_Descriptor<S> & d, int rc_ind = -1);
      void init_correct_map();

    public:

      S correct_value(const pcps::Configuration<S> & config, const char val_type) const;

      /// \brief function to access the pairing matrix data object by reference
      pcps::WfnData<S> & pairing_matrix_data() { return *this->_wf_deps->at(0); }

      /// \brief function to access the pairing matrix data object by constant reference
      const pcps::WfnData<S> & pairing_matrix_data() const { return *this->_wf_deps->at(0); }

      /// \brief function to access the pairing matrix by reference
      std::vector<S> & pairing_matrix() { return this->pairing_matrix_data().data_vec(); }

      /// \brief function to access the pairing matrix by constant reference
      const std::vector<S> & pairing_matrix() const { return this->pairing_matrix_data().data_vec(); }

      /// \brief function to access the sampling flag by reference
      bool & wfn_sampling() { return _wfn_sampling; }

      /// \brief function to access the sampling flag
      bool wfn_sampling() const { return _wfn_sampling; }

      /// \brief return the proposed ratio sum
      S prop_ratio_sum() const { return _prop_ratio_sum; }

      /// \brief return the value at the most recently initialized configuration
      pcps::MantExp<S> init_value() const { return _init_value; }

      /// \brief get the pqpq integrals from the Hamiltonian
      void get_pqpq_ints(const pcps::Hamiltonian<S> & ham);

      void pqpq_ratios(const pcps::Configuration<S> & config, std::vector<S> & ratios);

      RestPairingFunc(const pcps::Input & userinp,
                      pcps::Archive & arch,
                      const std::map<std::string, boost::shared_ptr<pcps::WfnData<int> > > & wi_map,
                      const std::map<std::string, boost::shared_ptr<pcps::WfnData<S> > > & wf_map,
                      const std::map<std::string, boost::shared_ptr<pcps::CoeffObj<S> > > & co_map);

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief converts the object into the command that will create the object
      ///
      /// \param[in,out]  arch     an archive in which to store the command
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void to_command(pcps::Archive & arch) const {
        arch << pcps::RestPairingFuncCommand::keyword() << std::endl;
        arch << this->_name                             << std::endl;
        arch << this->pairing_matrix_data().name()      << std::endl;
        arch << std::endl;
      }

      // pcps::MantExp<S> noiseless_value(const pcps::Configuration<S> & config);

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the value of the coefficient for an arbitrary configuration
      ///
      /// \param[in]      config   the configuration whose coefficient we want
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      pcps::MantExp<S> get_value(const pcps::Configuration<S> & config) {
        throw pcps::Exception("NOT YET IMPLEMENTED");
        return pcps::MantExp<S>();
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the value of the coefficient for the currently initialized configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      pcps::MantExp<S> get_value() {

        pcps::MantExp<S> retval;

        // if we are sampling with the wavefunction as the guiding function
        if ( _wfn_sampling ) {

          retval = this->_curr_val;

        // if we are sampling with the local rdm as the guiding function
        } else {

          retval = pcps::MantExp<S>(_acti_ratio_sum);

        }

        return retval;

      }

      pcps::MantExp<S> get_value(const pcps::Configuration<S> & new_config,
                                 const pcps::Configuration<S> & old_config,
                                 const pcps::ConfigDiff<S> & diff);

      void single_excitation_ratios(const pcps::Configuration<S> & config,
                                    std::vector<S> & a_ratios,
                                    std::vector<S> & b_ratios);

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Adds the derivatives of the coefficient with respect to the wavefunction's data
      ///        variables to the supplied vector after multiplying the derivatives by a constant.
      ///        The derivatives are computed assuming that the object is currently initialized for the
      ///        supplied lattice configuration.
      ///
      /// \param[in]      config  the lattice configuration the object is currently initialized for
      /// \param[in]      alpha   constant to multiply the coefficient's derivatives by
      /// \param[out]     der_vec vector to which to add the derivative contributions from this object
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void add_derivatives(const pcps::Configuration<S> & config,
                           const pcps::MantExp<S> & alpha,
                           std::vector<S> & der_vec)
      {
        throw pcps::Exception("NOT YET IMPLEMENTED");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief For each orbital pair excitation, contracts the wavefunction's derivatives at the
      ///        corresponding excited configuration with the supplied trial vector.
      ///        It is assumed that the object is initialized for the supplied configuration.
      ///
      /// \param[in]      config      the lattice configuration to excite from.
      /// \param[in]      alpha       constant to multiply the contraction by
      /// \param[in]      trial_vec   the vector to be contracted
      /// \param[in,out]  out_vec_a   vector to which to add the alpha orbital result
      /// \param[in,out]  out_vec_b   vector to which to add the beta orbital result
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void contract_der_to_orb_pairs(const pcps::Configuration<S> & config,
                                     const pcps::MantExp<S> & alpha,
                                     const std::vector<S> & trial_vec,
                                     std::vector<S> & out_vec_a,
                                     std::vector<S> & out_vec_b)
      {
        throw pcps::Exception("NOT YET IMPLEMENTED");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief For each wavefunction variable, contracts the various singly excited configurations'
      ///        wavefunction derivatives with respect to that variable with the supplied trial vector.
      ///        It is assumed that the object is initialized for the supplied configuration.
      ///
      /// \param[in]      config        the lattice configuration to excite from.
      /// \param[in]      alpha         constant to multiply the contraction by
      /// \param[in]      trial_vec_a   the vector to be contracted for alpha orbitals
      /// \param[in]      trial_vec_b   the vector to be contracted for beta orbitals
      /// \param[in,out]  out_vec       vector to which to add the result of the contraction
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void contract_orb_pairs_to_der(const pcps::Configuration<S> & config,
                                     const pcps::MantExp<S> & alpha,
                                     const std::vector<S> & trial_vec_a,
                                     const std::vector<S> & trial_vec_b,
                                     std::vector<S> & out_vec)
      {
        throw pcps::Exception("NOT YET IMPLEMENTED");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Multiply the supplied trial vector by the wavefunction derivative outer product matrix
      ///        (the local contribution to the overlap matrix) for the given configuration.
      ///        It is assumed that the object is initialized for the supplied configuration.
      ///
      /// \param[in]      config      the lattice configuration to excite from.
      /// \param[in]      alpha       constant to multiply the contraction by
      /// \param[in]      trial_vec   the vector to be contracted
      /// \param[in,out]  out_vec     vector to which to add the result of the contraction
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void contract_der_to_der(const pcps::Configuration<S> & config,
                               const pcps::MantExp<S> & alpha,
                               const std::vector<S> & trial_vec,
                               std::vector<S> & out_vec)
      {
        throw pcps::Exception("NOT YET IMPLEMENTED");
      }

      void initialize(const pcps::Configuration<S> & config);

      void update(const pcps::Configuration<S> & new_config,
                  const pcps::Configuration<S> & old_config,
                  const pcps::ConfigDiff<S> & diff);

  };

//  ///////////////////////////////////////////////////////////////////////////////////////////////////
//  /// \brief   Class to represent a restricted pairing function.  Here restricted means that each
//  ///          pair of spins must have an alpha and a beta, not that the overall wavefunction is
//  ///          necessarily a singlet.  If the pairing matrix is symmetric, however, then the
//  ///          the wavefunction is indeed a good singlet.
//  ///
//  ///          The pairing matrix rows correspond the alpha particles, while the columns correspond
//  ///          to beta particles.
//  ///
//  ///////////////////////////////////////////////////////////////////////////////////////////////////
//  template<class S> class RestPairingFunc : public pcps::CoeffObj<S> {
//
//    private:
//
//      /// \brief  the number of alpha sites
//      int _nas;
//
//      /// \brief  the number of alpha particles
//      int _nap;
//
//      /// \brief  the list of occupied alpha sites ordered to match the rows of the occupied pairing matrix
//      std::vector<int> _a_sites;
//
//      /// \brief  the list of occupied beta sites ordered to match the columns of the occupied pairing matrix
//      std::vector<int> _b_sites;
//
//      /// \brief  the row numbers of the occupied pairing matrix, ordered by their corresponding sites' indices
//      std::vector<int> _a_ordered_perm;
//
//      /// \brief  the column numbers of the occupied pairing matrix, ordered by their corresponding sites' indices
//      std::vector<int> _b_ordered_perm;
//
//      /// \brief  for each row of the occupied pairing matrix, the corresponding position in _a_ordered_perm
//      std::vector<int> _a_perm_pos;
//
//      /// \brief  for each column of the occupied pairing matrix, the corresponding position in _b_ordered_perm
//      std::vector<int> _b_perm_pos;
//
//      /// \brief  integer workspace vector
//      std::vector<int> _iwork;
//
//      /// \brief  the inverse of the occupied pairing matrix
//      std::vector<S> _inv_mat;
//
//      /// \brief  a matrix used as workspace
//      std::vector<S> _work_mat;
//
//      /// \brief  product of the inverse of the occupied pairing matrix with the 1st noise column vector
//      std::vector<S> _Iu1;
//
//      /// \brief  product of the inverse of the occupied pairing matrix with the 2nd noise column vector
//      std::vector<S> _Iu2;
//
//      /// \brief  product of the inverse of the occupied pairing matrix with the 1st noise row vector
//      std::vector<S> _v1I;
//
//      /// \brief  product of the inverse of the occupied pairing matrix with the 2nd noise row vector
//      std::vector<S> _v2I;
//
//      /// \brief  the 3x3 matrix involved in the matrix determinant lemma and the Woodbury matrix identity
//      std::vector<S> _little_mat;
//
//      /// \brief  workspace for diagonalizing the little matrix and performing some matrix multiplication
//      std::vector<S> _more_work;
//
//      /// \brief  Matrix to hold the inverse times U matrix involved in the Woodbury matrix identity.
//      ///         Note that what we actually store is the transpose of this matrix.
//      std::vector<S> _bigUI;
//
//      /// \brief  Matrix to hold the V times inverse matrix involved in the Woodbury matrix identity.
//      std::vector<S> _bigVI;
//
//      /// \brief  workspace for representing lattice and noise configurations
//      pcps::Configuration<S> _work_config;
//
//      /// \brief  workspace for representing differences between lattice and noise configurations
//      pcps::ConfigDiff<S> _work_diff;
//
//      void initialize_permutation(const std::vector<int> & sites, std::vector<int> & ordered, std::vector<int> & perm_pos);
//      void populate_site_lists(const pcps::Configuration<S> & config);
//      void build_occupied_matrix(const pcps::Configuration<S> & config, std::vector<S> & mat);
//      void compute_noise_products(const pcps::Configuration<S> & config);
//      int get_row_col_index(int site) const;
//      void build_little_mat(const pcps::ConfigDiff<S> & diff);
//      pcps::MantExp<S> get_single_excite_value(int new_ind, int old_ind, const std::vector<S> & m) const;
//      void get_row_col_diff(const pcps::ConfigDiff<S> & diff, S * const vec);
//      void print_occupancy_info();
//      //void swap_row_col(const int b, const int i, const int j);
//      void check_initialize(const pcps::Configuration<S> & config, const bool check);
//      void test_noise_removal(const pcps::Configuration<S> & config);
//      void get_single_excitation_a_ratios(const pcps::Configuration<S> & config, const int b, std::vector<S> & a_ratios);
//
//    public:
//
//      /// \brief function to access the pairing matrix data object by reference
//      pcps::WfnData<S> & pairing_matrix_data() { return *this->_wf_deps->at(0); }
//
//      /// \brief function to access the pairing matrix data object by constant reference
//      const pcps::WfnData<S> & pairing_matrix_data() const { return *this->_wf_deps->at(0); }
//
//      /// \brief function to access the pairing matrix by reference
//      std::vector<S> & pairing_matrix() { return this->pairing_matrix_data().data_vec(); }
//
//      /// \brief function to access the pairing matrix by constant reference
//      const std::vector<S> & pairing_matrix() const { return this->pairing_matrix_data().data_vec(); }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief Constructs the pairing function by reading it from an archive.
//      ///
//      /// \param[in]      userinp  the user's input options
//      /// \param[in,out]  arch     an archive from which to read the correlator
//      /// \param[in]      wi_map   map of names to available integer wavefunction data objects
//      /// \param[in]      wf_map   map of names to available floating point wavefunction data objects
//      /// \param[in]      co_map   map of names to available coefficient objects
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      RestPairingFunc(const pcps::Input & userinp,
//                      pcps::Archive & arch,
//                      const std::map<std::string, boost::shared_ptr<pcps::WfnData<int> > > & wi_map,
//                      const std::map<std::string, boost::shared_ptr<pcps::WfnData<S> > > & wf_map,
//                      const std::map<std::string, boost::shared_ptr<pcps::CoeffObj<S> > > & co_map)
//        : CoeffObj<S>(1, 0, 0),
//          _nas(userinp.nsites()/2),
//          _nap(userinp.nparticles()/2),
//          _a_sites(userinp.nparticles()/2, 0),
//          _b_sites(userinp.nparticles()/2, 0),
//          _a_ordered_perm(userinp.nparticles()/2, 0),
//          _b_ordered_perm(userinp.nparticles()/2, 0),
//          _a_perm_pos(userinp.nparticles()/2, 0),
//          _b_perm_pos(userinp.nparticles()/2, 0),
//          _iwork( 2 * std::max(3, userinp.nparticles()/2), 0),
//          _inv_mat(userinp.nparticles() * userinp.nparticles() / 4, pcps::zero<S>()),
//          _work_mat(userinp.nparticles() * userinp.nparticles() / 4, pcps::zero<S>()),
//          _Iu1(userinp.nparticles()/2, pcps::zero<S>()),
//          _Iu2(userinp.nparticles()/2, pcps::zero<S>()),
//          _v1I(userinp.nparticles()/2, pcps::zero<S>()),
//          _v2I(userinp.nparticles()/2, pcps::zero<S>()),
//          _little_mat(9, pcps::zero<S>()),
//          _more_work(std::max(9, 3*userinp.nparticles()/2), pcps::zero<S>()),
//          _bigUI(3*userinp.nparticles()/2, pcps::zero<S>()),
//          _bigVI(3*userinp.nparticles()/2, pcps::zero<S>()),
//          _work_config(userinp),
//          _work_diff(userinp)
//      {
//      
//        // check that we have zero net spin
//        if (userinp.z_spin() != 0.0)
//          throw pcps::Exception("the restricted pairing function is only defined for total Sz = 0");
//
//        // check that we an even particle number
//        if (userinp.nparticles() % 2 != 0)
//          throw pcps::Exception("the restricted pairing function is only defined for even particle number");
//
//        // get name of the pairing function
//        pcps::mpi::read_and_bcast(arch, this->_name, "Failed to read name in RestPairingFunc constructor.");
//
//        // get values
//        std::string wf_name;
//        pcps::mpi::read_and_bcast(arch, wf_name, "Failed to read name of value data in RestPairingFunc constructor.");
//        if (wf_map.count(wf_name) == 0)
//          throw pcps::Exception("no floating point data by the name %s") % wf_name;
//        this->_wf_deps->at(0) = wf_map.find(wf_name)->second;
//
//        // check that the value data object has the correct number of elements
//        if ( this->pairing_matrix().size() != userinp.nsites() * userinp.nsites() / 4 )
//          throw pcps::Exception("Wrong number of values in RestPairingFunc constructor.  Found %i values, expected %i.")
//                % this->pairing_matrix().size() % ( userinp.nsites() * userinp.nsites() / 4 );
//
//      }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief converts the object into the command that will create the object
//      ///
//      /// \param[in,out]  arch     an archive in which to store the command
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      void to_command(pcps::Archive & arch) const {
//        arch << pcps::RestPairingFuncCommand::keyword() << std::endl;
//        arch << this->_name                             << std::endl;
//        arch << this->pairing_matrix_data().name()      << std::endl;
//        arch << std::endl;
//      }
//
//      pcps::MantExp<S> noiseless_value(const pcps::Configuration<S> & config);
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief returns the value of the coefficient for an arbitrary configuration
//      ///
//      /// \param[in]      config   the configuration whose coefficient we want
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      pcps::MantExp<S> get_value(const pcps::Configuration<S> & config) {
//        throw pcps::Exception("NOT YET IMPLEMENTED");
//        return pcps::MantExp<S>();
//      }
//
//      pcps::MantExp<S> get_value(const pcps::Configuration<S> & new_config,
//                                 const pcps::Configuration<S> & old_config,
//                                 const pcps::ConfigDiff<S> & diff);
//
//      void single_excitation_ratios(const pcps::Configuration<S> & config,
//                                    std::vector<S> & a_ratios,
//                                    std::vector<S> & b_ratios);
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief Adds the derivatives of the coefficient with respect to the wavefunction's data
//      ///        variables to the supplied vector after multiplying the derivatives by a constant.
//      ///        The derivatives are computed assuming that the object is currently initialized for the
//      ///        supplied lattice configuration.
//      ///
//      /// \param[in]      config  the lattice configuration the object is currently initialized for
//      /// \param[in]      alpha   constant to multiply the coefficient's derivatives by
//      /// \param[out]     der_vec vector to which to add the derivative contributions from this object
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      void add_derivatives(const pcps::Configuration<S> & config,
//                           const pcps::MantExp<S> & alpha,
//                           std::vector<S> & der_vec)
//      {
//        throw pcps::Exception("NOT YET IMPLEMENTED");
//      }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief For each orbital pair excitation, contracts the wavefunction's derivatives at the
//      ///        corresponding excited configuration with the supplied trial vector.
//      ///        It is assumed that the object is initialized for the supplied configuration.
//      ///
//      /// \param[in]      config      the lattice configuration to excite from.
//      /// \param[in]      alpha       constant to multiply the contraction by
//      /// \param[in]      trial_vec   the vector to be contracted
//      /// \param[in,out]  out_vec_a   vector to which to add the alpha orbital result
//      /// \param[in,out]  out_vec_b   vector to which to add the beta orbital result
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      void contract_der_to_orb_pairs(const pcps::Configuration<S> & config,
//                                     const pcps::MantExp<S> & alpha,
//                                     const std::vector<S> & trial_vec,
//                                     std::vector<S> & out_vec_a,
//                                     std::vector<S> & out_vec_b)
//      {
//        throw pcps::Exception("NOT YET IMPLEMENTED");
//      }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief For each wavefunction variable, contracts the various singly excited configurations'
//      ///        wavefunction derivatives with respect to that variable with the supplied trial vector.
//      ///        It is assumed that the object is initialized for the supplied configuration.
//      ///
//      /// \param[in]      config        the lattice configuration to excite from.
//      /// \param[in]      alpha         constant to multiply the contraction by
//      /// \param[in]      trial_vec_a   the vector to be contracted for alpha orbitals
//      /// \param[in]      trial_vec_b   the vector to be contracted for beta orbitals
//      /// \param[in,out]  out_vec       vector to which to add the result of the contraction
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      void contract_orb_pairs_to_der(const pcps::Configuration<S> & config,
//                                     const pcps::MantExp<S> & alpha,
//                                     const std::vector<S> & trial_vec_a,
//                                     const std::vector<S> & trial_vec_b,
//                                     std::vector<S> & out_vec)
//      {
//        throw pcps::Exception("NOT YET IMPLEMENTED");
//      }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief Multiply the supplied trial vector by the wavefunction derivative outer product matrix
//      ///        (the local contribution to the overlap matrix) for the given configuration.
//      ///        It is assumed that the object is initialized for the supplied configuration.
//      ///
//      /// \param[in]      config      the lattice configuration to excite from.
//      /// \param[in]      alpha       constant to multiply the contraction by
//      /// \param[in]      trial_vec   the vector to be contracted
//      /// \param[in,out]  out_vec     vector to which to add the result of the contraction
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      void contract_der_to_der(const pcps::Configuration<S> & config,
//                               const pcps::MantExp<S> & alpha,
//                               const std::vector<S> & trial_vec,
//                               std::vector<S> & out_vec)
//      {
//        throw pcps::Exception("NOT YET IMPLEMENTED");
//      }
//
//      void initialize(const pcps::Configuration<S> & config);
//
//      void update(const pcps::Configuration<S> & new_config,
//                  const pcps::Configuration<S> & old_config,
//                  const pcps::ConfigDiff<S> & diff);
//
//  };

}

#endif
